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We demonstrate that the proper calculation of the linear response for finite-size systems can only 
be performed if the coupling to the leads/baths is explicitly taken into consideration. We exemplify 
this by obtaining a Kubo-type formula for heat transport in a finite-size system coupled to two 
thermal baths, kept at different temperatures. We show that the proper calculation results in a 
well-behaved response, without the singular contributions from degenerate states encountered when 
Kubo formulae for infinite-size systems are inappropriately used for finite-size systems. 
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The Liouville-von Neumann equation of motion for the 
density matrix p{t) of a closed system is {fi — 1): 



dp{t) 
dt 



LHp{t) = -i[H,p], 



(1) 



where H is the Hamiltonian describing the evolution of 
the system. U H = Hq + V, where F is a static weak 
coupling to an external field, one can use perturbation 
theory to find a solution p (t) = pQ + Sp (t) near a state 
Pq of the unperturbed system, LhqPo = 0. Neglecting 
LySpit) in Eq. ([T]), we integrate the resulting equation 
for 5p {t) to obtain the t ^ oo, steady-state solution: 



5p^ dte^'^'-'^'L 



vPo, 



(2) 



where rj — )■ 0+. If po = -^e~^^" describes the unper- 
turbed system in equilibrium at temperature ksT = 
then this is known as the Kubo formula [H, 0] : 



Sp = 



V{-t),. 



(3) 



where Vit) = g^iHaty^-tHot ^ -yy^g ^-^e identity 

[V{-t),e-f^"»] = -ie-^^o /(f drV {-t - ir) to rewrite: 
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5p. 



dte''^^ / drpoVi-t-ir) 
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In terms of the eigenbasis Ho\n) = e„ |n), {m\V {t) \n) 
i (cm - £«) V^nTie*^^'"""")*, leading to: 
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Note that there is no contribution from states with 
— Cm, for which (to|F (t) |n) = 0. This also fol- 
lows directly from Eq. ([3]); if we write V — Vb + Vj_, 
where Vq — J2t =e Vmn\iTi){n\ commutes with Hq, 
then [V{~t),po] = [V±{-t), po]. The "diagonal" part Vq 
of V does not contribute to 6p, and consequently has no 
influence on the static response functions. 



The lack of contributions from eigenstates with e„ — 
trn is, however, puzzling, because the well-known Drude 
weight derived directly from the Kubo formula is Q: 



D = ^-l y 



(m|J|n)p 



(6) 



i.e. it has contributions only from these eigenstates. 

To understand the reason for this difference, con- 
sider the derivation of Eq. © from Eq. (U), for sim- 
plicity for a one-dimensional chain described by Hq = 
+h.c.^ + VoJ2i^i^i+i^ where ni = c|q, 
plus a static potential V = J2i ^I'^i induced by a homoge- 
neous electric field E = —W. From the continuity equa- 

tion, v{t) = EiVifMt) = -^EiViiJi+iit) ~ Mt)], 

where Ji = it (c^j^Q — cj^ci+i^ is the local current opera- 
tor. This can be changed to X^/I^-i'^il*) "^''^'(O] = 
yVJ2iJi{t) = -EJ{t), where J{t) is the total cur- 
rent operator. Using V{t) — ~EJ{t) in Eq. ([4]) gives 
Sp = E dte~^* Jq drpoJ {—t - ir). The dc conduc- 
tivity is then a = dte~^* dr (J {—t — it) J), where 
(O) = Tr[poO], from which Eq. © follows. 

The only questionable step in this derivation, and the 
one responsible for going from a result with no contri- 
butions from states with e„ = e,,,, to one with contribu- 
tions only from these states, is the change J^i ^iJi+iit) ~^ 
J2i Vi-iJi{t)- This is only justified for an infinite system 
(where boundary terms are negligible), or a system with 
periodic boundary conditions and an external electric po- 
tential with the same periodicity f4]. It is certainly not 
valid for a finite size system connected to external leads, 
which break this symmetry. In such cases, the use of 
formulae like Eq. (|6]) is simply not appropriate iS]. 

The solution, however, is not the use of Eqs. (0]) or 
(O, derived for a finite size, closed system. Instead, one 
needs to derive their analogue for an open system coupled 
to leads. The reason is that p = po + Sp of Eq. ([5]) does 
not describe a non-equilibrium stationary state (NESS), 
in which transport of charge or heat through the sys- 
tem is possible. Instead, p is a first-order approximation 
of the thermal equilibrium state pth = e~^^/Z of the 
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full Hamiltonian H = Hq + V . Normally (for example 
if invariance to time reversal symmetry is not broken), 
no currents are generated in a thermal equilibrium state 
and therefore no steady-state transport through the finite 
system can be described by this approach [q]. 

To prove the above statement relating p to pth, we 
take {'m\V\n) = if e„i — e„, since as already dis- 
cussed, the "diagonal" part Vq of ^ does not contribute 
to transport. Consider then the eigenbasis Ii\n) — en\n), 
to first order perturbation in V. Since (mlVln) — 
for all Cm = En, we can apply the first order perturba- 
tion theory for non-degenerate states to all the states, 
whether degenerate or not, to find e„ = e„ + 0{V'^) and 
1^) = \n)+Em,e^^e^ ^^^H + OiV^) . This immedi- 
ately leads to pth^= Ej^e-^'-\h){n\ = po + Sp+0{V^), 
where 6p is indeed given by Eq. ([SJ. 

Therefore, we conclude that in order to describe an 
NESS with current fiow, one needs to go beyond viewing 
the finite size system of interest as a closed system, and to 
explicitly consider its connection to leads. We note that 
the effects of boundary conditions have been considered 
previously, eg in Refs. 0, Si- In particular, Ref. [8] 
derived a Kubo-like formula that takes into consideration 
cross-boundary currents in a stochastic approach. 

In this Letter we present an alternative deterministic 
formulation that explicitly considers the effects of cou- 
pling to leads (for charge transport) or thermal baths 
(for heat transport) on the state of the system. It is 
based on the Redfield equation [9] which describes the 
evolution of the projected density matrix for the central 
system of interest, obtained if we start from the Liouville- 
von Neumann equation for the total density matrix de- 
scribing the system-f leads/baths and use the projection 
technique [l3| to trace out the leads/baths. 

For simplicity, we assume coupling to thermal baths 
kept at temperatures T^/j^ = T± and investigate the 
heat transport in the resulting steady state. If AT -C T, 
this leads to a Kubo-like formula which replaces Eq. Q . 
This approach can be generalized straightforwardly to 
derive a Kubo-like formula for charge transport. 

The Redfield equation has the general form 



dpjt) 
dt 



= [LH + LL{TL)+LR{Tn)]p{t), 



(7) 



where Lnp — —i[H,p\, just like for an isolated system, 
while / jj. are additional terms that describe the effects 
of the left/right thermal baths (assumed to be in equilib- 
rium at their corresponding temperatures Tj^/b) on the 
evolution of the system. The expressions for ii/i? de- 
pend on the Hamiltonian H of the system and on its 
coupling to the baths (an example is provided below). 

If AT <^ T, we can Taylor expand ii/ij and re-arrange 
the Redfield equation to read: 



where Lb{T) = Lji{T) + Ll(T) is the contribution from 
the thermal baths if both are kept at the same temper- 
ature, while Lp(AT) collects the terms proportional to 
AT. Here we assume that H = Hq, i.e. that the thermal 
coupling does not induce an interaction V in the system. 
For charge transport such a term appears, and its Liou- 
villian Ly should be grouped together with Lp. 

Again, we are interested in the t — ^ oo, stationary state 
solution p of the above equation, i.e. 



Lp^O 



(9) 



which we assume to be unique for any value of AT. This 
means that L has a single zero eigenvalue, and all its 
other (transient) eigenvalues have a negative real part. 
Note that L|at=o = Lhq + Lb{T) has this property. In 
fact, one can show that in this case the i — ?• oo solution 
converges to the expected thermal equilibrium for the 
system held at temperature T, po — ^e^^^° ^10] . 

Eq. can be solved numerically to find this eigen- 
state. We call this solution p^x- However, our goal is to 
obtain a Kubo-like formula. This can be done by anal- 
ogy with the calculation for the closed system discussed 
in the beginning of this work. The first step is to sep- 
arate the Liouvillian L of Eq. ([S]) into a "large" plus a 
"small" part. There are two possible choices: either take 
the "large" part to be Tq^"* = L//„ +Lb{T) with the per- 
turbation Ai(i) = Lp{AT), or take t[,^^ = Lhq and let 
AL(2) =Lb{T) + Lp{AT) be the perturbation. 

We begin with the first choice. Tq^-* = Lho+Lb{T) has 
eigenvalues |Tq^^| and left /right eigenvectors {|£p)}, 
{|7?.p)}. As discussed, the unique steady-state solution 
of L^a^po = is po = ^e"^^". The deviation due 
to the perturbation AL^^^^ is obtained like in Eq. 



poo 

5p'i^=Y. rfte4>-.*|7^,)(/:,|AT(l: 



(1) 



Po- 



rn 



Note that the only divergent term, due to Lq^q = 0, dis- 
appears because (rolAL^^'po = (po|AT(i)po = 0. To 
see why, we start from Eq. ([3]), L{po + 5p) = 0, project 
it on (pol and keep terms only to first order, to find 

= (pol (l^'^ + AT(i)) (po + M = (pol AT(i)po since 
-^o^Vo = 0. As a result, Eq. (ITU)l has only regular contri- 
butions. We denote Po + Sp^^ — p^''. A similar approach 
using the eigenvalues and eigenvectors of Lhq + Lb{T) 
has been suggested in Ref. [1^], but for the Lindblad 
equation [14| instead of the Redfield equation. 

However, Eq. pUj) is difficult to use in practice; find- 



dpjt) 
dt 



[Lho + Lb{T) + Lp{AT)] p{t) - Lp{t), (8) j^g aU eigenstates of l'^q'' is a hard task unless the system 
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has an extremely small Hilbert space. A computationally 
simpler solution is obtained if we combine the eigenequa- 
tion Lp = with the constraint Trp = 1 into a regular 
system of coupled equations Lp — v, where, in matrix 
terms, L is defined by replacing the first row of the equa- 
tion Lp = by Tr (p) — 1, so that is a vector whose 
first element is 1, all remaining ones being 0. As a re- 
sult det (Z) 7^ while det (L) = 0. If solved numerically, 
Lp — V produces the expected solution p^x ■ 

We can also solve it to obtain a Kubo-like formula by 
dividing L = L^^^-hAL(i). A gain, the overbar shows that 



in matrix terms, L^^^ is obtained from L)^' by replacing 
its first row with Trp — 1, while AZ*-^' is obtained from 
AL*^"'^^ by replacing its first row with zeros. Then: 



(1) 



-[ZW]-iaZ(i) 



po- 



(11) 



This is much more convenient because inverting the non- 
singular matrix Zq^'' is a much simpler task than finding 
all the eigenvalues and eigenvectors of ig^-*. We have 
verified that both schemes produce identical results. 

The second option is to take l'^^^ — Lh„ and AL^^^ = 
Lb{T) + Lp{AT). In this case, we can still choose the 

(2) 

stationary solution associated with Lq to be the ther- 
mal equilibrium state at T, pQ = ^e^^^° . However, this 
solution is no longer unique, in particular the thermal 
equilibrium state corresponding to any other tempera- 

(2) 

ture satisfies LJ, 'po = 0. Expanding the analogue of Eq. 
(|3]) in the eigenbasis of Hq, we find: 



n.m 



X (2) 

^Pk 



p^'. Note that unlike of Eq. 



We call Po 

([TU]) , this solution has divergent contributions from states 
with e„ — £,„. As such, it is analogous to the Kubo 
formulae for infinite systems. 

This similarity is not accidental. Kubo formulae for 
infinite systems always ignore the coupling to the leads, 
taking Lq — Lhq and assuming that po — e^^^°/Z, 
where the temperature is arbitrarily chosen. Moreover, 
the driving force leading to transport is not a term 
Lp(AT), since the leads are ignored, but rather the ad- 
dition of some potential V to Hq leading to AL = Ly- 
Using only such a 1^ is rather questionable even for charge 
transport, because of the previously mentioned problems 
with the periodic boundary conditions, which may be 
negligible for infinite size systems but are certainly not 
for finite-size systems. For heat transport, using a F to 
describe a variation in the gravitational potential [2] is 
rather contrived, besides having the same boundary con- 
ditions issues. Nevertheless, with these assumptions, Eq. 
()12p maps into the usual Kubo formula for an infinite 
system. Thus, we can think of Eq. (IT^ with Lp{AT) 
included in AL as being a finite-size analog of the usual 



infinite-size Kubo formulae. (As discussed, for a finite 
size system Lp(AT) cannot be entirely replaced by an 
Ly, if steady-state transport is to be established). 

To see which of these two solutions - the regular solu- 
tion p^'' or the singular solution p^^ which is analogous 
to formulae for infinite systems - is the proper one, we 
compare them against the exact numerical solution p^^ 
of Eq. ® in the limit AT < T. 

We do this for a chain of N spins ^ coupled by nearest- 
neighbour exchange and placed in a magnetic field: 



N 



. (2) .-r-^ (m|Ai(2)p(,|,^\ 



i=l 2=1 

while the heat baths are collections of bosonic modes, 

where a — R,L indexes the riglit/left-side baths. The 
system-baths coupling is chosen as: 

where iL — 1, ip = N , i.e. the left/right bath is coupled 
to the first/last spin and can induce its spin-flipping. 

Using the projector technique [ll|, IHl , the equation of 
evolution for p(t) = TrBPt{t), where pt is the total den- 
sity of states for the system-hbaths, is found to second- 
order perturbation theory in Vint, to be: 



dt 



-i[no, p{t)] - J2 (K>ap(<)] 



a=L,R 



where 



(13) 

s\ -Sq. Here, (•) refers to the element-wise 
product of two matrices, {n\a ■ h\m) — {n\a\m) {n\h\m) . 
The bath matrices are defined in terms of the eigen- 
states of the system's Hamiltonian 7^o|^) = ^n\n) as: 



(a) |2 



S„ = TT^ \m){n\ e {VLynn)na {^mn) Da {^mn) \V^ 



+ e (r!„^) (1 + lla (nnrn)) {^nm) 1^^"^ P 



where fJ^ 



-r2„„i and fcm„ is defined by 



^kmn.a ^mm i-C is a bath mode resonant with this 



transition. Furthermore, Q{x) is the Heaviside function, 

na{^) = [e'^°^ — l] is the Bose-Einstein equilibrium 
distribution for the bosonic modes of energy il at the 
bath temperature T^ — I //3a, and Da{^) is the bath's 
density of states. The product Da {^mn) l^k"^ P is the 
bath's spectral density function. For simplicity, we take 
it to be a constant independent of m and n. 

Eq. (|13p is thus a particular example of the general Eq. 
([7]), with the last two terms coming from the coupling 
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FIG. 1: (color online) Di (squares) and D2 (circles) of Eq. 
dH} (a) vs. 5 = AT/2T at a fixed A = O.f, and (b) vs. A at a 
fixed S = 0.1, for iV = 8, J = 0.1, = 1. The insets show the 
steady-state thermal current calculated with pex (triangles) 
and p^'^' (squares and circles). See text for more details. 

to the two thermal baths. Since the temperatures Tj^/r 
enter only in the Bose-Einstein occupation numbers, it 
is straightforward to expand them when T/^//? = T ± 
AT < T and so to identify LgiT) and Lp(AT). 
We characterize the distance between the exact numer- 
ical solution pex and the two possible Kubo solutions p)^ , 
« = 1, 2 by calculating the norm: 




For the proper solution, this difference should be small 
but finite due to higher-order perturbation terms. 

Results typical of those found in all the cases we inves- 
tigated are shown in Fig. [T]for N = 8,J = 0.1, = 1.0. 
In Fig. (Ha), we plot Di,2 vs. 6 = AT/2T for a fixed 
system-baths coupling A = 0.1, while in Fig. [Ijb) we 
show them vs. A, for ST = Q.l. In both cases, D2 (cir- 
cles, axis on the right) is very large. In fact, because 
of the singular contributions from e„ — e„i states, D2 is 
divergent, its magnitude being controlled by the cutoff 77 
used (77 = 10~^ here). In contrast, Di (squares, left axis) 
has a small value independent of 77 0, showing that 
the regular p^'' is the proper Kubo solution. The insets 
show the thermal current calculated with p^x, p^^^ and 
Pj^ (triangles, squares, respectively circles). J2 becomes 
independent of 77 as 77 0, however, unlike Ji it is quite 
different from the exact solution. This again confirms 
that p^'' is the proper Kubo solution. 

Fig. [IJb) also shows that the answer depends on the 
details of the coupling to the baths. This is not sur- 
prising for a finite-size system: the intrinsic conductance 
of the system is added to comparable "contact" contri- 
butions from the interfaces between the system and the 
baths, and experiments measure the total conductance. 
It follows that quantitative modeling of transport in fi- 



nite systems will require a careful consideration of the 
entire experimental set-up. 

To conclude, we make two claims regarding the proper 
Kubo formulae to be used for finite-size systems. The 
first is that the coupling to leads/baths simply cannot 
be ignored, as is usually done for infinite-size systems. 
Instead, it has to be considered explicitly if steady-state 
transport is to be established in the system. Secondly, 
we showed that the proper way to derive a Kubo- like for- 
mula in such cases leads to a well-behaved result. This 
is to be contrasted with the formulae typically used in 
literature, similar to those valid for infinite-size systems, 
and which have singular contributions from eigenstates 
with e„ = In infinite systems this is not a prob- 
lem because the spectrum is continuous and the final re- 
sult is still well-behaved. However, for finite-size systems 
the spectrum is discrete and the singularities cannot be 
avoided. This is clearly unphysical: a finite-size system 
cannot have singular response functions. Of course, in re- 
ality the "eigenstates" of the open system are no longer 
sharp, instead they acquire a finite lifetime due to tunnel- 
ing of charge/heat into and out of the leads/baths. The 
need to properly consider the effects of the leads/baths 
on the system is therefore unavoidable. 
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